
################################################################################
# Troop-providers’ Ideational Commitment to UN Peacekeeping and Effectiveness ##
# Main replication code                                                       ##
# Journal of International Interactions                                       ##
# Author: Burak Giray, bgiray@uh.edu                                          ##
# Software used: RStudio version 1.4.1717                                     ##
################################################################################

library(readr)
library(plm)
library(stargazer)
library(ggplot2)
library(AER)
library(MASS)
library(sjPlot)
library(sjmisc)
library(tidyverse)
library(fixest)
library(simPH)
library(readr)
library("survival")
library("survminer")
library(readr)
library(stargazer)

########################
#### IMPORTING DATA ####
########################

duration_mydata <- read.csv("/Users/burakgiray/Desktop/Journal Submissions/Journal of International Interaction/replication files/duration_mydata.csv")
mydata <- read.csv("/Users/burakgiray/Desktop/Journal Submissions/Journal of International Interaction/replication files/mydata.csv")
PKOivc <- read_csv("/Users/burakgiray/Desktop/Journal Submissions/Journal of International Interaction/replication files/PKOivc.csv")

mydata$newviolence = mydata$civilianvio/100
mydata$lognewviolence = log(mydata$civilianvio + 1)
mydata$loglagnewviolence = log(mydata$lag_civilianvio + 1)
mydata$newtotalviolence = mydata$totalviolence/100
mydata = mydata %>% group_by(Mission) %>% arrange(date) %>% 
  mutate(lag_pko_fatalities = lag(pko_fatalities))

mydata = mydata %>% group_by(Mission_Country) %>% arrange(year) %>% 
  mutate(lag_polity2 = lag(polity2))

mydata = mydata %>% group_by(Mission) %>% arrange(date) %>% 
  mutate(lag_logtroopsize = lag(logtroopsize), lag_total_frac = lag(total_frac))
mydata2 =mydata %>% filter(!Mission_Country == "Rwanda")


##############################################
#### MODELS AND FIGURES IN THE MANUSCRIPT ####
##############################################

fit1 <- coxph(Surv(tstart, tstop, ep_end) ~ tot_composition+total_frac+logtroopsize +logtroopquality+logtroopfiscal+logDAC+pko_fatalities+logruralpop+ troopproviders+no_pa+polity2+logcivilianvio,data = duration_mydata)
summary(fit1)

test.ph = cox.zph(fit1) 
test.ph 

BaseVars <- c('logtroopfiscal')
duration_data1 <- tvc(duration_mydata, b = BaseVars, tvar = 'tstop', tfun = 'log')

# Table 1 of the Manuscript
solution1 <- coxph(Surv(tstart, tstop, ep_end) ~ tot_composition+total_frac+logtroopsize +logtroopquality+logtroopfiscal+logtroopfiscal_log+logDAC+pko_fatalities+ logruralpop+ troopproviders+no_pa+polity2+logcivilianvio,data = duration_data1)
summary(solution1)

## Table 2 of the Manuscript: 
main_one = femlm(civilianvio ~ tot_composition + total_frac+logtroopsize +logtroopquality+ logtroopfiscal+ logDAC+ pko_fatalities+ ruralpop+ troopproviders+difftime+ no_pa+polity2| Mission, mydata2, family = "negbin")
main_one_lagged = femlm(civilianvio ~ lag_civilianvio+tot_composition + total_frac+logtroopsize +logtroopquality+ logtroopfiscal+ logDAC+ pko_fatalities+ ruralpop+ troopproviders+difftime+ no_pa+polity2| Mission, mydata2, family = "negbin")
interaction_one = femlm(civilianvio ~ tot_composition + tot_composition:logtroopsize + total_frac+logtroopsize +logtroopquality+ logtroopfiscal+ logDAC+ pko_fatalities+ ruralpop+ troopproviders+difftime+ no_pa+polity2| Mission, mydata2, family = "negbin")
interaction_lagged = femlm(civilianvio ~ lag_civilianvio+tot_composition+ tot_composition:logtroopsize + total_frac+logtroopsize +logtroopquality+ logtroopfiscal+ logDAC+ pko_fatalities+ ruralpop+ troopproviders+difftime+ no_pa+polity2| Mission, mydata2, family = "negbin")

m.one = summary(main_one, se = "hetero") 
m.two = summary(main_one_lagged, se = "hetero") 
m.three = summary(interaction_one, se = "hetero") 
m.four = summary(interaction_lagged, se = "hetero") 

stargazer(coeftest(m.one),coeftest(m.two),coeftest(m.three),coeftest(m.four), no.space = TRUE, type="text")

# Figure 1 of the Manuscript
plot_interaction_one = summary(interaction_one, se = "hetero") 
summary(interaction_lagged, se = "hetero") 

library(MASS)
library(scales)
library(marginaleffects)
options(scipen=999)
plot = plot_cme(plot_interaction_one, effect = "tot_composition", condition = "logtroopsize", type = "link")+
  xlab("\n Troop Size (logged)") +
  ylab("Marginal Effect of Commitment to UNPKO\n")+
  scale_y_continuous(labels = label_number(scale = 1e+06/1000))+
  geom_hline(yintercept = 0, linetype = "dashed") +theme_minimal()
plot

###################################################
#### MODELS AND FIGURES IN THE ONLINE APPENDIX ####
###################################################

# Figure A1:
PKOivc1= PKOivc %>% group_by(Contributor) %>% summarize(ivc = mean(ivc))

PKOivc_highest = PKOivc1 %>% slice_max(ivc, n = 20)
ggplot(PKOivc_highest) +
  geom_point(aes(x = ivc, y = reorder(Contributor, ivc)))+	
  theme_grey() +xlab("") +ylab("")

PKOivc1 = PKOivc1 %>% filter(ivc != 0)
PKOivc_lowest = PKOivc1 %>% slice_min(ivc, n = 20)
PKOivc_lowest$Contributor[PKOivc_lowest$Contributor == "The former Yugoslav Republic of Macedonia"] = "North Macedonia"
PKOivc_lowest$Contributor[PKOivc_lowest$Contributor == "Central African Republic"] = "Central Afr. Rep."

ggplot(PKOivc_lowest) +geom_point(aes(x = ivc, y = reorder(Contributor, ivc)))+	theme_grey() +xlab("") +ylab("")

# Figure A2:
commitment_s = mydata %>% group_by(Mission) %>% summarize(commitment = mean(tot_composition))

scores_highest = commitment_s %>% slice_max(commitment, n = 20)
ggplot(scores_highest) +geom_point(aes(x = commitment, y = reorder(Mission, commitment)))+theme_grey() +xlab("") +ylab("")

scores_lowest = commitment_s %>% slice_min(commitment, n = 20)
ggplot(scores_lowest) +geom_point(aes(x = commitment, y = reorder(Mission, commitment)))+	theme_grey() +xlab("") +ylab("")

# Figure A3:
sample1 = duration_mydata %>% group_by(Mission) %>%  mutate(months_total = max(tstart), mean_comp = mean(tot_composition))
sample2 = sample1 %>% select(Mission, months_total, mean_comp) %>%  group_by(Mission) %>% distinct()
sample2$binary_comp = ifelse(sample2$mean_comp >= 0.55,1,0)
ended = duration_mydata %>% group_by(Mission) %>%  filter(ep_end == 1) %>% select(Mission, ep_end) %>% distinct()

sample3 <- merge(sample2,ended,by.x  = c("Mission"),by.y  = c("Mission"),all.x = TRUE,all.y =FALSE) %>%
  mutate_at(
    vars(ep_end), 
    ~ replace(
      ., 
      is.na(.), 
      0
    )
  )

kaplan_meier <- survfit(Surv(months_total, ep_end) ~ binary_comp,data = sample3)
ggsurvplot(kaplan_meier,conf.int = FALSE,risk.table.col = "strata", ggtheme = theme_bw(), palette = c("black", "red"),legend.labs =c("Low Commitment to UNPKO", "High Commitment to UNPKO"), xlim = c(0, 170))

# Table A3:
test.ph = cox.zph(fit1) 
test.ph 

# Figure A4:
set.seed(2)
Sim1 <- coxsimLinear(solution1, b = "tot_composition", Xj = seq(0.1, 1, by = 0.5),qi = "Relative Hazard",
                     nsim = 5, means = FALSE, ci = 0.90, extremesDrop = TRUE)
simGG(Sim1) +xlab("Commitment to UNPKO")

# Table A4:
fit2 <- coxph(Surv(tstart, tstop, ep_end) ~ tot_composition+total_polarization+logtroopsize +logtroopquality+logtroopfiscal+logDAC+pko_fatalities+ logruralpop+ troopproviders+no_pa+polity2+logcivilianvio,data = duration_mydata)
test.ph2 = cox.zph(fit2)
test.ph2
solution2 <- coxph(Surv(tstart, tstop, ep_end) ~ tot_composition+ total_polarization+logtroopsize +logtroopquality+logtroopfiscal+logtroopfiscal_log+logDAC+pko_fatalities+ logruralpop+ troopproviders+no_pa+polity2+logcivilianvio,data = duration_data1)

size_added_frac <- coxph(Surv(tstart, tstop, ep_end) ~ tot_composition + total_frac+logtroopsize+logpolicesize+logciviliansize +logtroopquality+ logtroopfiscal+ logDAC+ pko_fatalities+ logruralpop+ troopproviders+ no_pa+polity2+ logcivilianvio, data = duration_mydata)
test.ph3 = cox.zph(size_added_frac)
test.ph3

size_added_pol <- coxph(Surv(tstart, tstop, ep_end) ~ tot_composition + total_polarization+logtroopsize+logpolicesize+logciviliansize +logtroopquality+ logtroopfiscal+ logDAC+ pko_fatalities+ logruralpop+ troopproviders+ no_pa+polity2+ logcivilianvio,data = duration_mydata)
test.ph4 = cox.zph(size_added_pol)
test.ph4

full <- coxph(Surv(tstart, tstop, ep_end) ~ tot_composition + total_frac+logtroopsize+logpolicesize+logciviliansize +logtroopquality+ logtroopfiscal+ logDAC +lognonDAC+logEU+ pko_fatalities+ logruralpop+ troopproviders+ no_pa+polity2+ logcivilianvio,data = duration_mydata)
test.ph5 = cox.zph(full)
test.ph5

full2 <- coxph(Surv(tstart, tstop, ep_end) ~ tot_composition + total_polarization+logtroopsize+logpolicesize+logciviliansize +logtroopquality+ logtroopfiscal+ logDAC +lognonDAC+logEU+ pko_fatalities+ logruralpop+ troopproviders+ no_pa+polity2+ logcivilianvio,data = duration_mydata)
test.ph6 = cox.zph(full2)
test.ph6
BaseVars <- c('logEU')
duration_data4 <- tvc(duration_mydata, b = BaseVars, tvar = 'tstop', tfun = 'log')
full2solution <- coxph(Surv(tstart, tstop, ep_end) ~ tot_composition + total_polarization+logtroopsize+logpolicesize+logciviliansize +logtroopquality+ logtroopfiscal+ logDAC +lognonDAC+logEU +logEU_log+ pko_fatalities+ logruralpop+ troopproviders+ no_pa+polity2+ logcivilianvio,data = duration_data4)

stargazer(solution2, size_added_frac, size_added_pol,full,full2solution, no.space = TRUE, type="text") 


# Figure A5:
heteroscedasticity = glm.nb(civilianvio ~ tot_composition+logtroopsize +logtroopquality+ logtroopfiscal+ logDAC+ pko_fatalities+ ruralpop+ troopproviders+difftime+ no_pa+polity2, mydata2)
summary(heteroscedasticity)

par(mfrow = c(2, 2))
plot(heteroscedasticity)
lmtest::bptest(heteroscedasticity)

# Figure A6:
install.packages("corrplot")
library(corrplot)
my_data = mydata2 %>% select(tot_composition, total_frac, total_polarization, logtroopquality, logtroopfiscal) %>% ungroup()
my_data = my_data %>% select(-Mission)
my_data  = my_data %>% rename(Polarization = total_polarization, Fractionalization = total_frac, `Commitment to UNPKO` = tot_composition, `Troop quality`=logtroopquality, `Fiscal capacity`= logtroopfiscal) 
my_data <- my_data[complete.cases(my_data), ]

res <- cor(my_data)
testRes = cor.mtest(res, conf.level = 0.95)
corrplot(res, p.mat = testRes$p, method = 'color', diag = FALSE, type = 'upper',
         sig.level = c(0.001, 0.01, 0.05), pch.cex = 0.9,
         insig = 'label_sig', pch.col = 'grey20', order = 'AOE')


# Table A5:
bivariate_robustse = femlm(civilianvio ~ tot_composition | Mission, mydata2, family = "negbin") #noRwanda
bivariate_robustse2 = femlm(civilianvio ~ tot_composition | Mission, mydata, family = "negbin") # withRwanda
br_one = summary(bivariate_robustse, se = "hetero")
br_two = summary(bivariate_robustse2, se = "hetero")

stargazer(coeftest(br_one),coeftest(br_two), no.space = TRUE, type="text")

# Figure A7:
plot(mydata2$tot_composition, mydata2$civilianvio, pch = 19, col = "black")

ggplot(data.frame(mydata2),aes(x=tot_composition,y=civilianvio))+
  geom_point() + 
  geom_smooth(method = "glm.nb") + xlab("Commitment to UNPKO") + ylab("Civilian Victimization")


# Table A6:
no_independentvariable1 = femlm(civilianvio ~  total_frac+logtroopsize +logtroopquality+ logtroopfiscal+ logDAC+ pko_fatalities+ ruralpop+ troopproviders+difftime+ no_pa+polity2| Mission, mydata2, family = "negbin")
no_independentvariable_frac = femlm(civilianvio ~  total_frac+logtroopsize +logtroopquality+ logtroopfiscal+ logDAC+ pko_fatalities+ ruralpop+ troopproviders+difftime+ no_pa+polity2| Mission, mydata, family = "negbin")
no_independentvariable2 = femlm(civilianvio ~  total_polarization+logtroopsize +logtroopquality+ logtroopfiscal+ logDAC+ pko_fatalities+ ruralpop+ troopproviders+difftime+ no_pa+polity2| Mission, mydata2, family = "negbin")
no_independentvariable_pol = femlm(civilianvio ~  total_polarization+logtroopsize +logtroopquality+ logtroopfiscal+ logDAC+ pko_fatalities+ ruralpop+ troopproviders+difftime+ no_pa+polity2| Mission, mydata, family = "negbin")

no_in1 = summary(no_independentvariable1, se = "hetero") # no rwanda
no_in2 =summary(no_independentvariable_frac, se = "hetero") 
no_in3 =summary(no_independentvariable2, se = "hetero") # no rwanda
no_in4 =summary(no_independentvariable_pol, se = "hetero") 

stargazer(coeftest(no_in1),coeftest(no_in2),coeftest(no_in3),coeftest(no_in4), no.space = TRUE, type="text")

# Table A7:
b3 = femlm(civilianvio ~ tot_composition + total_frac+logtroopsize +logtroopquality+ logtroopfiscal+ logDAC+ pko_fatalities+ ruralpop+ troopproviders+difftime+ no_pa+polity2| Mission, mydata, family = "negbin")
b3_lagged = femlm(civilianvio ~ lag_civilianvio+tot_composition + total_frac+logtroopsize +logtroopquality+ logtroopfiscal+ logDAC+ pko_fatalities+ ruralpop+ troopproviders+difftime+ no_pa+polity2| Mission, mydata, family = "negbin")
b3_interaction = femlm(civilianvio ~ tot_composition + tot_composition:logtroopsize + total_frac+logtroopsize +logtroopquality+ logtroopfiscal+ logDAC+ pko_fatalities+ ruralpop+ troopproviders+difftime+ no_pa+polity2| Mission, mydata, family = "negbin")
b3_interaction_lagged = femlm(civilianvio ~ lag_civilianvio+tot_composition+ tot_composition:logtroopsize + total_frac+logtroopsize +logtroopquality+ logtroopfiscal+ logDAC+ pko_fatalities+ ruralpop+ troopproviders+difftime+ no_pa+polity2| Mission, mydata, family = "negbin")

fulls1 = summary(b3, se = "hetero") 
fulls2 =summary(b3_lagged, se = "hetero") 
fulls3 =summary(b3_interaction, se = "hetero") 
fulls4 =summary(b3_interaction_lagged, se = "hetero") 

stargazer(coeftest(fulls1),coeftest(fulls2),coeftest(fulls3),coeftest(fulls4), no.space = TRUE, type="text")

# Table A8:
b4 = femlm(civilianvio ~ tot_composition + total_polarization+logtroopsize +logtroopquality+ logtroopfiscal+ logDAC+ pko_fatalities+ ruralpop+ troopproviders+difftime+ no_pa+polity2| Mission, mydata2, family = "negbin")
b4_lagged = femlm(civilianvio ~ lag_civilianvio+tot_composition + total_polarization+logtroopsize +logtroopquality+ logtroopfiscal+ logDAC+ pko_fatalities+ ruralpop+ troopproviders+difftime+ no_pa+polity2| Mission, mydata2, family = "negbin")
b4_interaction = femlm(civilianvio ~ tot_composition + tot_composition:logtroopsize + total_polarization+logtroopsize +logtroopquality+ logtroopfiscal+ logDAC+ pko_fatalities+ ruralpop+ troopproviders+difftime+ no_pa+polity2| Mission, mydata2, family = "negbin")
b4_interaction_lagged = femlm(civilianvio ~ lag_civilianvio+tot_composition+ tot_composition:logtroopsize + total_polarization+logtroopsize +logtroopquality+ logtroopfiscal+ logDAC+ pko_fatalities+ ruralpop+ troopproviders+difftime+ no_pa+polity2| Mission, mydata2, family = "negbin")

pol1 = summary(b4, se = "hetero") 
pol2 = summary(b4_lagged, se = "hetero") 
pol3 = summary(b4_interaction, se = "hetero") 
pol4 = summary(b4_interaction_lagged, se = "hetero") 

stargazer(coeftest(pol1),coeftest(pol2),coeftest(pol3),coeftest(pol4), no.space = TRUE, type="text")


# Table A9:
b5 = femlm(civilianvio ~ tot_composition + total_frac+logtroopsize+logpolicesize+logciviliansize +logtroopquality+ logtroopfiscal+ logDAC +lognonDAC+logEU+ pko_fatalities+ ruralpop+ troopproviders+difftime+ no_pa+polity2| Mission, mydata2, family = "negbin")
b5_lagged = femlm(civilianvio ~ lag_civilianvio+tot_composition + total_frac+logtroopsize+logpolicesize+logciviliansize +logtroopquality+ logtroopfiscal+ logDAC +lognonDAC+logEU+ pko_fatalities+ ruralpop+ troopproviders+difftime+ no_pa+polity2| Mission, mydata2, family = "negbin")
b5_interaction = femlm(civilianvio ~ tot_composition + tot_composition:logtroopsize + total_frac+logtroopsize+logpolicesize+logciviliansize +logtroopquality+ logtroopfiscal+ logDAC +lognonDAC+logEU+ pko_fatalities+ ruralpop+ troopproviders+difftime+ no_pa+polity2| Mission, mydata2, family = "negbin")
b5_interaction_lagged = femlm(civilianvio ~ lag_civilianvio+tot_composition+ tot_composition:logtroopsize + total_frac+logtroopsize+logpolicesize+logciviliansize +logtroopquality+ logtroopfiscal+ logDAC +lognonDAC+logEU+ pko_fatalities+ ruralpop+ troopproviders+difftime+ no_pa+polity2| Mission, mydata2, family = "negbin")

controls1 = summary(b5, se = "hetero") 
controls2 = summary(b5_lagged, se = "hetero") 
controls3 = summary(b5_interaction, se = "hetero") 
controls4 = summary(b5_interaction_lagged, se = "hetero") 

stargazer(coeftest(controls1),coeftest(controls2),coeftest(controls3),coeftest(controls4), no.space = TRUE, type="text")

# Table A10:
# Autocracy variable
mydata2$autocracy = NULL
mydata2$autocracy[mydata2$polity2 <= -6 ] = 1
mydata2$autocracy[is.na(mydata2$autocracy)] = 0 
mydata2$autocracy[is.na(mydata2$polity2)] =NA

## Autocracy table : 

autoc_one = femlm(civilianvio ~ tot_composition + total_frac+logtroopsize +logtroopquality+ logtroopfiscal+ logDAC+ pko_fatalities+ ruralpop+ troopproviders+difftime+ no_pa+autocracy| Mission, mydata2, family = "negbin")
autoc_one_lagged = femlm(civilianvio ~ lag_civilianvio+tot_composition + total_frac+logtroopsize +logtroopquality+ logtroopfiscal+ logDAC+ pko_fatalities+ ruralpop+ troopproviders+difftime+ no_pa+autocracy| Mission, mydata2, family = "negbin")
autoc_interaction_one = femlm(civilianvio ~ tot_composition + tot_composition:logtroopsize + total_frac+logtroopsize +logtroopquality+ logtroopfiscal+ logDAC+ pko_fatalities+ ruralpop+ troopproviders+difftime+ no_pa+autocracy| Mission, mydata2, family = "negbin")
autoc_interaction_lagged = femlm(civilianvio ~ lag_civilianvio+tot_composition+ tot_composition:logtroopsize + total_frac+logtroopsize +logtroopquality+ logtroopfiscal+ logDAC+ pko_fatalities+ ruralpop+ troopproviders+difftime+ no_pa+autocracy| Mission, mydata2, family = "negbin")

a.one = summary(autoc_one, se = "hetero") 
a.two = summary(autoc_one_lagged, se = "hetero") 
a.three = summary(autoc_interaction_one, se = "hetero") 
a.four = summary(autoc_interaction_lagged, se = "hetero")

stargazer(coeftest(a.one),coeftest(a.two),coeftest(a.three),coeftest(a.four), no.space = TRUE, type="text")


## Table A11: 
logged1 = plm(lognewviolence ~ tot_composition + total_frac+logtroopsize +logtroopquality+ logtroopfiscal+ logDAC+ pko_fatalities+ ruralpop+ troopproviders+difftime+ no_pa+polity2, mydata2, index = c("Mission"), model = "within")
logged2 = plm(lognewviolence ~ loglagnewviolence+tot_composition + total_frac+logtroopsize +logtroopquality+ logtroopfiscal+ logDAC+ pko_fatalities+ ruralpop+ troopproviders+difftime+ no_pa+polity2, mydata2,  index = c("Mission"), model = "within")
logged3 = plm(lognewviolence ~ tot_composition + tot_composition:logtroopsize + total_frac+logtroopsize +logtroopquality+ logtroopfiscal+ logDAC+ pko_fatalities+ ruralpop+ troopproviders+difftime+ no_pa+polity2, mydata2,  index = c("Mission"), model = "within")
logged4 = plm(lognewviolence ~ loglagnewviolence+tot_composition+ tot_composition:logtroopsize + total_frac+logtroopsize +logtroopquality+ logtroopfiscal+ logDAC+ pko_fatalities+ ruralpop+ troopproviders+difftime+ no_pa+polity2, mydata2,  index = c("Mission"), model = "within")

stargazer(logged1,logged2,logged3,logged4, no.space = TRUE, type="text")


## Table A12:
robustness_nested <- plm(tot_composition ~ loglagnewviolence +lag_pko_fatalities+ laglogDACT1  +lagruralT1+laggdpT1+ troopproviders+no_pa+difftime+polity2, 
                         data = mydata,
                         index = c("Mission"), 
                         model = "within")

robustness_nested2 <- plm(tot_composition ~ loglagnewviolence +lag_pko_fatalities+ laglogDACT1  +lagruralT1+laggdpT1+ troopproviders+no_pa+difftime+polity2, 
                          data = mydata2,
                          index = c("Mission"), 
                          model = "within")

r.one =coeftest(robustness_nested2, vcov. = vcovHC, type = "HC1")
r.two =coeftest(robustness_nested, vcov. = vcovHC, type = "HC1")

stargazer(r.one,r.two, no.space = TRUE, type="text")














